Assessment of interatomic potentials for atomistic analysis 
of static and dynamic properties of screw dislocations in W 



D. Cereceda and J. M. Perlado 

Universidad Politecnica de Madrid 
28005 Madrid, Spain 

A. Stukowski, S. Queyreau, and J. Marian* 

Lawrence Livermore National Laboratory 
Livermore, CA 94551 

Lisa Ventelon and M.-C. Marinica 

Service de Recherches de Metallurgique Physique, CEA-DEN, F-91 191, Gif-sur-Yvette, 
France 

M. R. Gilbert 

EURATOM/CCFE Fusion Association, Culham Science Centre, Abingdon, OX14 3DB, 
United Kingdom 

E-mail: *marianl@llnl . gov 

Abstract. Screw dislocations in bcc metals display non-planar cores at zero temperature 
which result in high lattice friction and thermally activated strain rate behavior. In bcc W, 
electronic structure molecular statics calculations reveal a compact, non-degenerate core with 
an associated Peierls stress between 1.7 and 2.8 GPa. However, a full picture of the dynamic 
behavior of dislocations can only be gained by using more efficient atomistic simulations based 
on semiempirical interatomic potentials. In this paper we assess the suitability of five different 
potentials in terms of static properties relevant to screw dislocations in pure W. As well, we 
perform molecular dynamics simulations of stress-assisted glide using all five potentials to 
study the dynamic behavior of screw dislocations under shear stress. Dislocations are seen to 
display thermally-activated motion in most of the applied stress range, with a gradual transition 
to a viscous damping regime at high stresses. We find that one potential predicts a core 
transformation from compact to dissociated at finite temperature that affects the energetics 
of kink-pair production and impacts the mechanism of motion. We conclude that a modified 
embedded-atom potential achieves the best compromise in terms of static and dynamic screw 
dislocation properties, although at an expense of about ten-fold compared to central potentials. 
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1. Introduction. 

Tungsten (W) is being considered as a leading candidate for plasma-facing applications in 
magnetic fusion energy (MFE) devices. The most attractive properties of W for MFE are its 
high melting point and thermal conductivity, low sputtering yield and low long-term disposal 
radioactive footprint. These advantages are accompanied unfortunately with very low fracture 
toughness characterized by brittle trans- and inter- granular failure, which severely restrict the 
useful operating temperature window [Q]. 

Transgranular plasticity in refractory metals, including W, is governed by the temperature 
dependence of screw dislocation motion. W is typically alloyed with 5~26 at.% Re to increase 
low temperature ductility and improve high temperature strength and plasticity J2]|. The 
physical origins behind the Re-induced ductilization have been discussed in the literature 
OHO and point in some way or another to alterations in the core structure of | (111) screw 
dislocations, which both reduce the effective Peierls stress o> and extend the number of 
possible slip pathways. A direct consequence of a reduced Peierls stress, e.g. as via Re 
alloying, is an enhanced dislocation mobility at low temperatures. Recent electronic structure 
calculations of a P in pure W give values between 1.7 and 2.8 GPa [|3l |6]|. This means that, 
under most conditions relevant to technological applications, where stresses are of the order 
of only a few hundred MPa, a reduction in op of a few hundred MPa may not be significant 
to the plastic behavior of W and W alloys. Instead, it is the thermally-activated and three- 
dimensional character of screw dislocation motion, the associated solution softening behavior, 
as well as the temperature dependence of the core structure, that control bulk ductility. 

All of these aspects cannot be studied in atomistic detail using current experimental 
capabilities. By contrast, atomistic methods based on semiempirical potentials have enabled 
large-scale molecular dynamics (MD) simulations, so that, at present, calculations of single- 
dislocation mobility, core structure and transformations, etc., can be obtained with reasonable 
accuracy. However, care must be exercised when choosing from the dozen or so W potentials 
available in the literature. Semiempirical force fields with both pair and cohesive contributions 
(e.g. following the embedded atom method formalism) are typically considered to achieve an 
optimum balance between efficiency and accuracy. These are typically fitted to reproduce 
some basic bulk and defect properties such as lattice parameter, elastic constants, vacancy 
formation energy, surface energies, etc., but generally not dislocation properties. Of these, 
it is known that the screw dislocation core structure at K should be non-degenerate (also 
known as compact), as revealed by density functional theory (DFT) calculations 01 El 

Previous atomistic calculations on screw dislocations in W have been performed by 
Mrovec et al [7], Fikar et al [8] and Tian and Woo [9J. Mrovec et al studied the dislocation 
core structure and calculated the Peierls stress at K using a tight-binding-based bond-order 
potential (TB-BOpfl. They predicted a non-degenerate core structure and a Peierls stress of 
4.3 GPa. For their part, Fikar et al studied core structures and energies of screw dislocations 
using three different interatomic potentials, all of which display dissociated cores. Lastly, Tian 

* BOP potentials include non-central atomic interactions to represent the effect of (^-electrons in transition 
metals 
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and Woo examined the mobility of screw dislocations also with an embedded-atom potential 
that predicts a dissociated core structure. They were able to obtain dislocation velocities at 
stresses above the Peierls stress at K. However, no systematic study of dislocation motion 
in W at finite temperature has been conducted. Characterizing dislocation motion in the 
stress-temperature space is important to parameterize the so-called mobility functions used 
in higher-level methods such as dislocation dynamics. The purpose of the mobility functions 
is to provide a quantitative measure of the response of dislocations to applied and internal 
stresses. 

Unfortunately, one of the most important difficulties associated with such studies is the 
scale-dependent nature of MD simulations, which require exceedingly high strain rates to 
drive the system over time scales accessible computationally, of the order of a a few tens 
of ns. Because of these limitations, MD is incapable at present of properly capturing the 
thermally-activated motion of screw dislocations at low stresses. However, MD simulations 
can still provide valuable input in intermediate-to-high stress conditions and in situations 
where the deformation rates are high. The objective of this paper is to compare five different 
interatomic potentials -that have not been fitted against screw dislocation data- and assess 
their performance in terms of static and dynamic screw dislocation properties. By static 
properties we mean several reference parameters at K as obtained with DFT calculations. 
The dynamic behavior is evaluated in terms of screw dislocation mobility as a function of 
stress and temperature. Due to the absence of 'reference' mobility data against which to 
compare the potentials, we will simply draw several general conclusions based on the inter- 
comparison among potentials. 

The paper is organized as follows. First, we discuss the distinctive features of each 
potential and calculate the structure of a screw dislocation core. The Peierls potential and 
the 7 surface are then calculated and verified against existing DFT and TB-BOP calculations. 
Subsequently, we introduce the computational setup for the dynamic mobility simulations and 
calculate dislocation velocities as a function of temperature and stress. Subsequently, a study 
of the core trajectories in the plane defined by the glide and normal directions is carried out. 
We finish by analyzing the causes of the temperature-dependent behavior of each potential and 
emphasizing the insufficiency of static calculations to fully characterize dislocation motion at 
finite temperatures. 

2. Computational details 

2.1. Interatomic potentials 

Our calculations have been performed with the parallel MD code LAMMPS U5J. Table 
CD gives basic information about the five different potentials considered here, among which 
there are three embedded-atom method (EAM) potentials, one Tersoff-Brenner-type bond- 
order potential (TF-BOP), and one modified EAM (MEAM). Note that the TB-BOP used by 
Mrovec et al was deemed not suitable for dynamics simulations by its authors [10J and 
has thus not been considered here. Hereafter they are referred to in the text by the identifiers 
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given in the table header. This selection of W potentials, from the dozen or so available in the 
literature, is not meant to be an implicit assessment of the quality of those not employed here. 

Table 1. Properties of potentials used: lattice parameter ao, shear modulus \i, Peierls stress 
up, computational cost, core structure at K, and thermal expansion coefficient a. Potentials 
EAM1, EAM2 and TF-BOP display a threefold symmetric (degenerate) core, while EAM3 
and MEAM predict compact (non-degenerate) cores. The values of the volumetric thermal 
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ma. 

Two important quantities for characterizing screw dislocation cores at K are the Peierls 
potential, defined as the energy path from one equilibrium position to another on a {110} 
plane, and the 7 surface along the [111] direction also on {110} planes. The Peierls potential 
governs the morphology of kinks (e.g. [1221 ) while workers such as Duesbery and Vitek ||23l 
have provided evidence for a direct correspondence between the shape of the |(111){110} 
gamma surface and the screw dislocation core structure. These are plotted, respectively, for 
all potentials in Figs. [Hand [2] on the (110) plane. DFT data for both calculations are also 
shown for comparison. 

The Peierls potential was obtained using the nudged elastic band method [1241 in the 
manner proposed by Groger and Vitek [|25l . whereas the DFT calculations in both cases were 
obtained using a plane wave self-consistent field code as described in Refs. [|26l and 11271 . 
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Figure 1. Peierls potential for all potentials tested here. DFT calculations from Ventelon et al 
ll27l are shown for comparison. 
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Figure 2. (110) 7 surface along the [111] direction for all the potentials considered in this 
work. DFT calculations from Ventelon et al ll27l are shown for comparison. 



2.2. Simulation setup 

To measure dislocation velocities, we have performed stress-controlled simulations of 
7; (111) dislocations with the maximum resolved shear stress (MRSS) on a {112} plane. The 
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justification to focus on {112}-type planes is twofold. First, as Argon and Maloof [fTTTl have 
shown, under tensile loading, most loading orientations and temperatures, result in some 
degree of |(111){112} slip. Also, Li et al [|5]| have shown that {112} slip is important in 
W alloys with high Re concentrations. Second, certain EAM potentials intrinsically deviate 
from MRSS behavior when the MRSS plane is of the {110} type for reasons that have been 
discussed at length in the literature [Q2] l3Tl[i~3l . 

We have provided the details regarding the computational setup for this type of 
simulations in prior publications |[T3l[T4| . Suffice it to say that orthogonal boxes of sufficient 
size with axes x = ^ [111], y = [110], and z = [112], corresponding to the line, glide, and 
plane normal directions, respectively, were used to mitigate finite size effects (cf. Ref. IfTBl ). 
a xz is then applied on the computational cell boundaries and simulations are conducted in the 
NPT ensemble. Periodic boundary conditions are used in the line and glide directions. The 
reference cell dimensions were L x = 25 ^a , L y = 100 [v^ao] > and L z = 50 [\/6ao] , 
where the amounts in brackets are the dimensions of the nominal unit cell in the coordinate 
system employed here. The reference configuration contains 7.5 x 10 5 atoms, which results 
in strain rates of 1.4 x 10 6 ~ 7 s -1 for dislocation velocities between 10 and 100 m-s -1 . 

All simulations were run on LLNL's ATLAS cluster using 128 and 256 processors at a 
reference cost of ~ 10~ 7 CPU seconds per atom per time step for potential EAM1. 



3. Results 



The simulation setup, boundary conditions, and velocity calculations, as they relate to the 
present work, are discussed in depth by Cereceda et al lfi~4~l . The temperature and stress ranges 
covered were, respectively, 300 to 2100 K and 200 to 2000 MPa. All the simulations were run 
for 100 ps and configuration data were extracted every picosecond. The procedure to extract 
dislocation velocities from MD simulations is well established in the literature [|28l |29l [i~3l : 
from the position of the core, velocities are calculated as the derivative of the displacement- 
time curves for each case. 



3.1. Screw dislocation mobility 

Figure |3] shows all the (a-v ) data for the five interatomic potentials tested. The figures 
also contain the temperature dependence for each case. Generally, the velocities increase 
monotonically with stress and temperature, although at different rates depending on the 
potential. To first order, the mechanism of motion followed by the dislocations depends on the 
Peierls stress. This means that, at a maximum applied shear stress of 2000 MPa, the EAM1, 
EAM2, and MEAM potentials both operate under o> (cf. Tabled]), while for the EAM3 and 
TF-BOP there are several data points above it. In either case, dislocation motion is mostly 
governed by the thermally activated kink-pair nucleation mechanism, and thus display an 
exponential dependence with a and T. This can be qualitatively appreciated in the figure, 
although in the Appendix a more quantitative analysis is carried out. 

Another important aspect of dislocation motion is the extent of MRSS motion displayed, 
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Figure 3. Screw dislocation velocity as a function of applied shear stress and temperature for 
the five interatomic potentials considered here. Note that the velocity and stress axes are not 
on the same scale for each case. 



i.e. whether there are deviations from glide on the MRSS {112}-type plane. In Figure |4] 
we analyze the trajectories on the yz-plane for different combinations of a and T over 100 
ps of simulation. Perfect MRSS behavior is characterized by trajectories parallel to 0°. As 
the figure shows, all the EAM potentials display nearly perfect MRSS behavior, while for 
the MEAM small deviations in the acceleration phase are captured. The TF-BOP potential 
displays the most erratic motion with an overall deviation of the order of five degrees. At 
higher temperatures and stresses, this effect is enhanced to the point that the dislocation exits 
the simulation box only a few picoseconds after the shear stress is applied. This is the reason 
why there are fewer data points in the a-v curves shown in Fig. |3]for the TF-BOP potential. 
The trajectories shown in the figure are effective, i.e. they are not sufficiently time resolved to 
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capture the atomistic details of dislocation motion. Nevertheless, the operating mechanism of 
motion is by way of nucleation and propagation of kink pairs on {110} planes adjacent to the 
MRSS (112) plane. 
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Figure 4. Line-averaged dislocation trajectories on the yz-plane for two combinations of a 
and T and over 100 ps. Planes forming 0, 5, 15 and 30° with the (112) MRSS plane are 
represented with dotted lines (angles not to scale). Except for the TF-BOP potential, all the 
simulations yield small <5° deviations from MRSS motion. 





3.2. Dislocation core structure at finite temperature 

As shown in Fig.[3l the a-v data are not conducive to comparison among potentials. Instead, 
in Fig. [5] they are plotted as a function of interatomic potential for a number of selected 
temperatures. The figure reveals an interesting trend: the relative behavior of all the potentials 
remains unchanged for all temperatures with the exception of EAM3. At low temperatures, 
this potential exhibits a relatively high dislocation mobility, akin to that displayed by 'fast' 
potentials such as EAM2. However, above 900 K, the mobility is reduced (relative to the 



Cereceda et al 



9 



other interatomic potentials) to values more in line with 'slower' potentials such as EAM1. 
Moreover, if one examines the trajectories followed by the dislocation at 500 MPa0, a notable 





Figure 5. Comparison of interatomic potentials for the data given in Fig. [3] The colored 
lines correspond to exponential fits obtained in the Appendix for potentials EAM1, EAM3 and 
MEAM. 



difference in behavior within the EAM3 potential can be observed. At a temperature of 600 K, 
the dislocation follows a biased path on an effective glide plane forming f^30° with the MRSS 
plane. However, at 1800 K, the dislocation follows a path that deviates only slightly from that 
dictated by the Peach-Kohler force (i.e 0°). This is quantitatively displayed in Figured where 
this time the trajectories are resolved with atomistic detail. The figure shows unequivocally 
that dislocation motion proceeds via the formation of kink pairs on { 1 10 } planes bordering the 
MRSS [112] plane (at ±30°). Moreover, the details of the trajectory at 600 K suggest biased 
formation on the (101) plane (+30°), whereas at 1800 K random-walk behavior is displayed, 
with kink pairs forming equally on both available {110} planes. 

The behaviors illustrated in Figs. |5]and|6]for_potential EAM3 suggest a change in core 
structure with temperature for a given stress state Q To examine the physical structure of the 
dislocation core at different temperatures one can use time-averaged differential displacement 
(DD) maps (these maps were used in Table [T]for each K configurations). The DD maps are 
obtained by running MD simulations of crystals containing four screw dislocations arranged 

" To be meaningful, this analysis must be performed at relatively low stresses to interfere the least amount 
possible with the investigated temperature effect. 

+ We know that stress also induces its own core transformations as explained in Ref. lETTl . 
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EAM3 potential, a=500 MPa 




y-displacement [nm] 

Figure 6. Line-averaged dislocation trajectories on the yz-plane for the EAM3 potential at 500 
MPa. Results for two 200-ps temperatures above and below the presumed core transformation 
temperature of around 1200 K are shown. Planes forming 0, 5, 15 and 30° with the (112) 
MRSS plane are represented with dotted lines (angles not to scale). 

in a balanced quadrupole configuration and periodic boundary conditions. The size of the 
simulation box is 20 x 15 x 18 multiples of the bcc lattice vectors [111] x [121] x [101]. 
The dimensions are adjusted to the equilibrium lattice constant at the given temperature. 
For the finite temperature simulations, the displacement of each atomic string is determined 
by averaging over all 40 atoms in the string and over a time window of 100 fs, being 
sufficiently long to avoid noise due to thermal vibrations yet short enough to not capture 
diffusive behavior. The results are shown in Fig. [7]for configurations in the < T < 2100 
K interval. The figure confirms that the the EAM3 core is the only one showing an 
appreciable transformation from non-degenerate to degenerate, clearly seen at and above 
1500 K. Although DD maps are a useful tool to quickly analyze core structures, next we 
complement the results in Fig. [7] with a more quantitative approach based on fundamental 
lattice properties. 

3.3. Analysis of screw dislocation core stability. 

Duesbery and Vitek [23J have provided a simple rule that relates the shape of the | (111) { 110} 
7-surface to the core structure at K. They used the following inequality: 
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Figure 7. Time-averaged core structures for the five potentials tested here in the entire 
temperature range. 



to predict whether a screw dislocation will display a compact core. 7 (|) and 7 (|) are the 
energies corresponding to the | and | magnitudes of the generalized fault vectors, which can 
be obtained by reference to Fig. |2] The idea is that, if the above inequality is satisfied, |-type 
faults will be preferred over | ones, leading to non-dissociated core structures. However, 
although Duesbery and Vitek applied this simple rule to six different bcc metalsQ with 

* V, Cr, Nb, Mo, Ta, and W, all described by Finnis-Sinclair potentials 02 
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remarkable success, we find that in our case it does not hold for potentials EAM1 and EAM2. 
Thus, here we try a different approach based on the analysis carried out by Gilbert and 
Dudarev [|32l . 

These authors have shown that, alternatively, the screw dislocation core structure 

in bcc systems can be related to the periodic interaction energy between adjacent (111) 
strings of atoms forming the crystal. Their analysis, which was performed primarily to 
help guide potential development, provides a framework to predict whether the favored 
core structure at K is compact or dissociated. In particular, they derive the so-called 
"first-nearest-neighbour (INN) inter-string interaction law" of the potential and use this in 
a 2D Frenkel-Kontorova (2D-FK) model of interacting (111) strings to find the minimum 
energy screw core-structure. Here we extend their methodology to finite temperatures using 
the quasiharmonic approximation, i.e. by relating volume changes to temperature via pre- 
computed thermal expansion coefficients for each potential. In this fashion, we first compute 
the inter-string interaction laws as a function of the lattice parameter and then obtain the 
equivalent temperature as: T = 3 (a/do — 1) /a. Here, a is the thermal expansion coefficient 
(given for each potential in Table [TJ), a the lattice parameter at K, and a the lattice 
parameter corresponding to a temperature T (within the quasiharmonic approximation). For 
the reminder of this section, we refer to a as a(T) to highlight this temperature dependence. 

For each temperature T the INN inter-string interaction law U\ (d) was derived by rigidly 
translating a single (111) string with respect to a perfect lattice with lattice parameter a(T) 
and measuring the associated variation in energy under the particular interatomic potential. 
The resulting curve, which, according to the 2D-FK model defined in Ref. 11321 . is dominated 
by the contributions from the moving string interacting with its six INNs, can be unfolded 
using a Fourier analysis to produce the required pair-wise interaction law for the 2D-FK 
model. An example of such a law for EAM3 at K is shown in figure (8^b). A perfect screw 
dislocation, inserted into a lattice of (111) atomic strings (for a given a = a(T)) according 
to the isotropic elasticity solution, was then relaxed using U\(d) and the nature of the relaxed 
core was determined by visual inspection of its differential displacement map. 

Figure Oa) shows the variation in the favored core structure as a function of T for each 
of the five potentials. On the y-axis of the plot we have calculated the ratio of the string 
separation d* associated with the inflection points in the corresponding U\ (d) law (highlighted 
by the vertical dotted lines in figure [Hb)) to 5(T)/6, where b(T) is the corresponding Burgers 
vector of each potential as a function of T. 

As observed by Gilbert and Dudarev fl32l . the favored core structure depends on the 
position of the inflection points of the Ui(d) function. Specifically, a fully compact core 
is characterized by minimum string separations of b(T)/6, which are the in-line separation 
distances between each of the three (111) strings immediately surrounding the core (red 
circles in Figure |8{c), which shows a differential displacement map for a compact core) and 
their two nearest strings forming the next shell of strings out from the core (blue circles in 
[8tc)). When the inflection points in U\(d) are located at a distance of less than or equal to 
b(T)/6, then the compact, non-dissociated core is always stable. Furthermore, even if the 
separations d* associated with U"(d) = are such that \d*\ is somewhat greater than b(T)/6, 
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the compact core may still be stable provided that Z7i is only slowly varying around these 
inflection points — meaning that the forces (—U[(d)) between strings are relatively constant 
over a range of d values. This, for example, is the situation in the case of the MEAM potential 
in Fig. (Hfa), where the ratio d* / (6(T)/6) is greater than one for all T, but the favored core 
remains compact. 

However, when the ratio is significantly greater than one, as is the case for both EAM1 
and EAM2 at all temperatures, then the compact core becomes unstable and the secondary 
strings out from the core tend to move towards (along (111)) one of their primary- string 
neighbors (signified by the major arrows in the "arms" of the non-compact core shown for 
EAM1 in Tabled]) — ultimately leading to the stabilization of the non-compact, three-fold 
symmetric dissociated core. 

Thus, if, as a function of T, there is significant variation in this ratio, then the preferred 
equilibrium core structure can also change. In our analysis, we find that, consistent with the 
transition observed in Fig. [7J there is a large shift for EAM3 in the value of d* / (6(T)/6) {d* is 
such that U"(d*) = 0) above ~ 1500 K. At this point, the equilibrium core structure diverges 
from the compact core, and becomes more and more dissociated as temperature increases 
further. In the next section we discuss the implications of our findings. 

4. Discussion 

4.1. Comparison of static properties of potentials 

The five interatomic potentials tested here follow different formulations and have been fitted 
to different physical properties. It is not our objective to discuss the fitting process or the 
quality of each one, but only to discuss their performance in relation to screw dislocation 
modeling. When selecting potentials for screw dislocation simulations, two of the properties 
most looked at are the core structure and the Peierls stress. For W, these have been obtained 
using electronic structure calculations of different sorts, which reveal a compact core and ap 
between 1.7 and 2.8 GPa. Regarding the core structure at K, only EAM3 and MEAM 
reproduce it correctly, although, as shown in Section 13.21 the EAM3 potential does not 
preserve this structure at high temperature. In terms of op, the five potentials studied here 
give a range of values from 1.1 to 4.0 GPa. As Table Q] shows, EAM2 and EAM3 display 
values of 4.0 and 1.8 GPa, respectively. These potentials are also the most computationally 
efficient, which is an advantage when computing resources are limited. It is worth noting 
that the cost of these potentials was evaluated using the cutoff radii specified in the original 
references shown in Table CD 

Furthermore, the subspace of the energy landscape most relevant to screw dislocation 
motion is the Peierls potential and the 7 surface (Figs. CD and |2]). The shape of the Peierls 
potential is only correctly reproduced by potentials EAM3 and MEAM, while the rest predict 
trajectories with metastable states along their path, in disagreement with DFT calculations. 
Regarding the 7 surface, potentials EAM1, EAM3 and MEAM all predict the essential 
qualitative and quantitative features of the DFT results and are also in good agreement with 
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Figure 8. (a) Analysis of the favored core structure for the five different potentials as a function 
of temperature. A compact core is designated by a circle, and a non-compact by a diamond. 
The y-position of each point is the magnitude of distance of the inflection points in the INN 
string-interaction law (i.e. d* , which is the distance at which U"(d) = 0) normalized to the 
quantity b(T)/6, where b(T) = >/3a(T) /2. The dashed horizontal line separates the region of 
phase space where compact cores are favored (below) from the region where non-compact are 
favored (above), (b) The INN string-interaction law for EAM3 at K. The red dashed vertical 
lines indicate the position of the inflection points in this curve, (c) Differential displacement 
map of the compact core predicted at K for the U\ function from EAM3. Each of the three 
strings closest to the core (red circles) are separated from their two closest secondary strings 
(blue circles) by 6/6. In the figure the temperature dependence of b is omitted for the sake of 
clarity. 



the results by Groger et al using a TB-BOP OTTl . 

Thus, on the basis of all these calculations, the MEAM potential appears to be the best 
suited of those tried here to carry out dislocation simulations at any temperature. When 
computational cost is of the essence, EAM3 may be considered an acceptable replacement 
for static calculations or at low temperatures and stresses. 

4.2. Mobility of screw dislocations 

Dislocation mobility is highly multidimensional in that it displays multiparametric 
dependencies, e.g. on stress, temperature, dislocation character, slip system, etc. Dislocation 
velocities are difficult to infer from straining experiments, while they are costly and subjected 
to size limitations in simulations. Measurements ll35l and calculations ll36l l37ll of edge 
dislocation velocities have been carried out in W. However, other than the values computed 
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by Tian and Woo [9| at very high stress (>3.6 GPa), to our knowledge no data exist on screw 
dislocation mobility in W at low stresses. In this work, we have focused on the temperature 
and stress dependences, while we have kept the dislocation character and the slip system fixed. 

A quick look at Fig. [5] reveals several interesting details. First, the EAM2 and MEAM 
consistently give the highest and lowest velocities, respectively, regardless of temperature. 
Since atomistic simulations commonly overestimate screw dislocation velocities, particularly 
in the low-stress range, this may be another reason in favor of using MEAM. This is likely 
to be due to the fact that the dislocation core remains compact in the entire temperature 
range using the MEAM potential (cf. Fig. [7]). Second, screw dislocations move by thermally- 
activated mechanisms below the Peierls stress, transitioning to a viscous damping regime 
above it. At the maximum applied stress of 2000 MPa, some dislocations have been driven 
past the Peierls stress as given by their respective potentials (cf. Table [I])- This is certainly 
the case for the TF-BOP and possibly potentials EAM2 and EAM3. One would therefore 
expect to see a gradual exhaustion of the thermally-activated regime and a transition into a 
linear regime. Interestingly, such a transition appears to occur for potential EAM1 which has 
ap = 4.0 GPa. It was shown by Gilbert et al lfi~3~1 . however, that the actual transition stress 
decreases with the square of the temperature, which may be what is seen here. Appropriate 
exponential fits to the data shown in Fig. [5] carried out in the Appendix reveal useful 
parameters that define the thermally activated regime. 

4.3. Dislocation core transformation 

The behavior that emanates from the results in Figs. |5] and |6]for potential EAM3 is the 
manifestation of a temperature-driven dislocation core transformation that occurs as a result 
of changes to the free energy landscape. We have characterized this transformation via 
differential displacement maps of time-averaged atomic positions at finite temperatures, and 
a quasiharmonic analysis of the location of inflection points in the (111) interaction energy, 
which is known to control the dislocation core structure (cf. Figs. [6] and [8]). 

This is seen to affect the dislocation mobility as well. We have shown that the reported 
core transformation has an impact on both the stress and temperature dependencies. Indeed, 
in the analysis carried out in the Appendix, it is shown that the temperature dependence of 
fitted a-v relations at T < 1200 K for EAM3 does not carry over to higher temperatures. By 
contrast, the same analysis does not yield significant differences between the low and high 
temperature regimes for potentials EAM1 and MEAM. This is further indication that the core 
structure may impact the motion mechanisms in the corresponding temperature range. 

We emphasize, however, that as long as there does not exist independent evidence of 
this dislocation core structural change with temperature, the discussion about its true impact 
on the dynamic behavior of screw dislocations remains solely speculative, and we cautiously 
warn against using the EAM3 potential above the observed transformation temperature of 
~1500 K. In this sense, the quasi-harmonic analysis performed in Section [3731 would be very 
amenable to DFT calculations, as it consists solely of zero-temperature calculations. This 
would provide an independent means to prove or disprove -at least within the limitations of 
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the quasi-harmonic analysis- the behavior predicted by EAM3 at high temperature. 
4.4. Mechanism of motion 

The dislocation trajectories shown in Fig. @] demonstrate that screw dislocations move 
primarily along the direction of the applied stress. The only notable exception is the TF- 
BOP, for which significant transitions out of plane are observed. The figure, however, does 
not provide insights into the atomistic mechanism of motion. Then, in Fig. [6] trajectories for 
the EAM3 were analyzed with higher spatial and temporal resolution at temperatures 600 and 
1800 K. At both temperatures, the dislocation moves by elementary {110} kink-pair episodes. 
It is reasonable to assume that this mechanism can be extrapolated to other potentials that 
yield similar effective trajectories (close to MRSS plane). However, for the EAM3 results, 
there are some differences in terms of the temperature at which the trajectory was extracted. 
At 1800 K it appears as though the unit mechanism is composed of one +30° jump (on a (011) 
plane) followed by a correlated —30° jump (on a (101) plane). In other words, the dislocation 
appears to move by kink-pair episodes on the (112) plane that consist of two alternating and 
correlated ±30° kinks. This is consistent with the mechanism proposed by Duesbury 11381 . 
Overall, this results in a trajectory that follows a random walk and that, on average, forms 
zero degrees with the MRSS plane. Interestingly, at 600 K kink pairs on the +30° plane seem 
to be favored in a proportion of three or more to one over —30° ones. It is unclear at this point 
if the dislocation core transition discussed above for EAM3 is responsible for this difference. 
Again, as stated in Section l43l we are reluctant to construe this as real physical behavior until 
more is known about the core structure transformation. Our main message from the analysis 
of trajectories is that despite the MRSS plane being of the {112} family, motion proceeds by 
way of kink pairs on {110} planes presumably for all potentials. 

5. Summary 

To summarize, the main findings of this paper can be condensed into the following main 
items: 

• We have calculated static properties relevant to screw dislocations using five different 
interatomic potentials for W. These include three EAM, one BOP and one MEAM. 

• We have calculated screw dislocation mobilities for all potentials on a {112} glide plane. 
Our calculations provide elements to judge the MEAM potential as the most suitable for 
dislocation calculations. 

• We have observed a temperature-induced dislocation core transformation -from compact 
to dissociated- for one of the potentials tested. Lacking independent confirmation, we 
cannot confirm whether this corresponds to a real physical phenomenon or is an artifact, 
but the transformation is indeed seen to impact the dynamic properties of dislocations. 

• Our analysis of the five interatomic potentials suggests, first, that the atomistic nature 
of the dislocation core governs behavior at larger scales and, second, a purely static 
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treatment of the dislocation core is insufficient to predict and describe the dynamics of 
dislocations. 
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Appendix A. 

Here we analyze the overall impact of the dislocation core transition for potential EAM3 on 
dislocation mobility. We fit the data given in Fig.[3]to the general expression: 

v(a,T) = Aaexp (- H °~J V * ^ (A.l) 

where A, H , and V* are fitting constants that represent, respectively, a velocity prefactor, 
the kink pair energy at K, and the activation volume. We obtain these for three potentials, 
namely, EAM1, EAM3, and MEAM, and carry out the fit first including all temperatures. 
The results for each case are shown in Table IA1I These values deserve some commentary, 

Table Al. Fitting parameters for the analytical mobility function IA.1I The average fitting 
error for A, Ho, and V* was, respectively, 6, 9 and 10%. Regular script: values from full 
temperature fits; bold script: values from low temperature (< 900 K) fits; in parentheses: 
percentage difference between both sets of fits. 



Potential 


A 


ms^MPa -1 




#o [eV] 


V* [b 3 ] 


EAM3 

MEAM 

EAM1 


0.26 0.24 (8%) 
0.19 0.17(10%) 
0.30 0.28 (7%) 


0.05 0.04 (20%) 
0.08 0.07 (12%) 
0.08 0.07 (12%) 


0.42 0.19 (55%) 
0.26 0.23(12%) 
0.24 0.24 (0%) 



particularly H and V*. Using the method described by Ventelon and Willaime [26], we 
have obtained a kink-pair energy of H =1.7 eV for the MEAM potential. Experimentally, 
Brunner [43] has obtained a value of 1.75 eV from the temperature dependence of flow stress 
measurements in W, in very good agreement with the calculated value but significantly higher 
than the MD values. Giannattasio et al 11441 have obtained values of the order of 1.0 eV 
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inferred from the strain rate dependence of the brittle-to- ductile transition, still much larger 
than those reported here. Similarly, Tarleton and Roberts [|45l have found values of V*=20b z 
to be representative of the kink-pair process in W. Again, these are two orders of magnitude 
larger than ours. The magnitudes of H and V* obtained in our analysis suggest a very 
'soft' thermally activated process, something not necessarily consistent with the static data 
presented in Section |2T1 The low values of H and V* obtained in our simulations are likely 
due to overdriven screw dislocation dynamics in the MD simulations. 

Next, we obtain additional fits using only data at 300, 600, and 900 K, i.e. at temperatures 
below the presumed core transformation for potentail EAM3. The corresponding parameter 
values are shown in bold script for each case in Table lAll The percentage difference between 
the values for full and low temperature fits is given in parentheses. Although the differences 
in the parameter A are similar in all cases, those for H and V* are clearly largest for EAM3. 
Particularly, the differences in the activation volume, which is known to be most sensitive to 
the core structure, are considerable. This further reinforces the notion that, from a dynamical 
standpoint, dislocations are behaving differently below and above ^1000 K. In contrast, for 
the other two potentials, the dynamic behavior above and below this presumed transition 
temperature is governed by the same laws. The low temperature fits obtained here are shown 
in Fig. [51 The fits provide very good agreement with the EAM1 and ME AM data at all 
temperatures, whereas they gradually worsen for EAM3 as temperature increases. As well, 
the deviation of the fits above op for potential EAM3 can be clearly appreciated. 
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